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We introduce a theoretical approach to determine the spin structure of harmonically trapped 
atoms with two-body zero-range interactions subject to an equal mixture of Rashba and Dresselhaus 
spin-orbit coupling created through Raman coupling of atomic hyperfine states. The spin structure 
of bosonic and fermionic two-particle systems with finite and infinite two-body interaction strength g 
is calculated. Taking advantage of the fact that the IV-boson and iV-fermion systems with infinitely 
large coupling strength g are analytically solvable for vanishing spin-orbit coupling strength k so 
and vanishing Raman coupling strength fl, we develop an effective spin model that is accurate to 
second-order in fl for any k so and infinite g. The three- and four-particle systems are considered 
explicitly. It is shown that the effective spin Hamiltonian, which contains a Heisenberg exchange 
term and an anisotropic Dzyaloshinskii-Moriya exchange term, describes the transitions that these 
systems undergo with the change of k so as a competition between independent spin dynamics and 
nearest-neighbor spin interactions. 

PACS numbers: 


I. INTRODUCTION 


Spin-orbit (more precisely, spin-momentum) coupled 
systems continue to attract a great deal of attention 
due to the rich physics of the spin-Hall effect, topo¬ 
logical insulators, and Majorana fermions p]-(E|. While 
these topics fall traditionally into condensed matter ter¬ 
ritory, recent advances in cold atomic gases have led to a 
fruitful cross fertilization of atomic and condensed mat¬ 
ter physics. On the experimental side, synthetic gauge 
fields have been realized in ultracold atom systems [6H15] . 
On the theoretical side, systems with different kinds of 
spin-orbit coupling have been studied at the many- and 


few-body levels [6L I16l - t24| . Understanding the effects of 
spin-orbit coupling in few-atom systems opens the door 
to a bottom-up understanding of many-body systems. 
Nowadays, ultracold few-atom systems can be prepared 
and manipulated in experiments [25j, [26| . For example, 
taking advantage of a Fano-Feshbach resonance, the in¬ 
teraction between the atoms can be tuned [27|. More¬ 
over, the confinement geometry can be changed from 
three-dimensional to effectively two-dimensional to ef¬ 
fectively one-dimensional [28i] . These experimental ad¬ 
vances were guided by and stimulated a good number 
of theoretical few-body studies [2§-[42j|. Much analyti¬ 
cal work has been done by approximating the true alkali 
atom-alkali atom potential by a zero-range contact po¬ 


tential [29|, 1321434 1361 , 140M45| | . This approximation cap¬ 
tures the low energy regime reliably but fails to repro¬ 
duce high-energy properties (such as the characteristics 
of deeply bound states). Assuming zero-range interac¬ 
tions, the energy spectra of two one-dimensional bosons 
and two one-dimensional fermions with arbitrary two- 
body coupling constant and spin-orbit coupling strength 
have been calculated and the interplay between the spin- 
orbit coupling term, Raman coupling term, and the two- 
body potential has been analyzed [2(| . 


Extending earlier work 0, |2lJ, this paper investi¬ 
gates the spin structure of harmonically trapped atoms 
with two-body zero-range interactions subject to an 
equal mixture of Rashba and Dresselhaus spin-orbit cou¬ 
pling created through Raman coupling of atomic hyper- 
fine states. The spin structure of two identical one- 
dimensional bosons and two identical one-dimensional 
fermions with finite and infinite two-body interaction 
strength g is calculated for weak to strong spin-orbit 
coupling strength. The bosonic and fermionic systems 
display, in general, different behaviors as the spin-orbit 
coupling strength k so and the two-body coupling con¬ 
stant g are changed. For infinite g , however, the bosonic 
and fermionic systems display the same spin structure. 
An interesting transition of the spin structure is found in 
going from small to large k so . To understand the system 
behaviors for finite and infinite g, an effective Hamilto¬ 
nian for any k so that is accurate up to second order in the 
Raman coupling strength fl is derived. Effective Hamil¬ 
tonian have been utilized in various areas of physics and 
the Hubbard model [46| , the Born-Oppenheimer approxi- 
matio n [f7l , the Ising model [48j , and the Heisenberg spin 
chain [49| are prominent examples. Various approaches 
to generate effective Hamiltonian and the connection be¬ 
tween perturbation theory and selected effective Hamil¬ 
tonian have been discussed in Refs. [5CK55| . 

For infinitely large g and arbitrary N , we integrate out 
the spatial degrees of freedom and recast the resulting 
effective Hamiltonian in terms of spin operators. Single 
spin terms are proportional to fl while spin-spin interac¬ 
tions are proportional to fl 2 . Effective spin Hamiltonian 
have been derived previously for one-dimensional systems 
without spin-orbit and Raman coupling [56l-l59|. In those 
cases, spin-spin interactions were introduced by allowing 
for small deviations from \g\ = oo; in essence, this in¬ 
troduces a tunneling term. In our case, the single spin 
and the spin-spin terms are introduced by the Raman 
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coupling. The single spin term is proportional to f 1 and 
has been discussed in Ref. [2l|. In essence, the Raman 
coupling creates an effective local R-field that the spins 
follow. The spin-spin interaction term has, to the best of 
our knowledge, not been discussed before in this context. 
This term contains two contributions. The first contri¬ 
bution is of the type of the “usual” Heisenberg exchange 
term [45[, i.e., it is a aj ■ Ok term, where dj is the spin 
operator of spin j. The second contribution is of the type 
of the anisotropic Dzyaloshinskii-Moriya exchange term, 
i.e., it is a one-dimensional analog of the D ■ (dj x ak) 
term where D is a constant vector. 

The remainder of this paper is organized as follows: 
Section El defines the system Hamiltonian, discusses 
its symmetries and introduces an effective low-energy 
Hamiltonian. Section Ell determines the spin structure 
of the two-particle system for different g and k so . Using 
the effective low-energy Hamiltonian, Sec. IIVI derives an 
effective spin Hamiltonian for infinitely large g and ap¬ 
plies this Hamiltonian to determine the spin correlations 
of the three- and four-particle systems. A transition from 
a regime where the dynamics is governed by single spin 
physics to a regime where the dynamics is governed by 
spin-spin interactions is obtained. Finally, Sec. [V] sum¬ 
marizes and concludes. The appendices contain a number 
of technical details, including the evaluation of the ma¬ 
trix elements for infinite g , which are needed to construct 
both the effective Hamiltonian and the full Hamiltonian. 
The analytical techniques and results are expected to be 
useful for other one-dimensional studies. 


II. SYSTEM HAMILTONIAN AND GENERAL 
CONSIDERATIONS 

We consider N one-dimensional atoms of mass m in 
a harmonic trap with angular trapping frequency w and 
zero-range two-body interactions g8{xj — Xk), where Xj is 
the position of the jth particle with respect to the center 
of the trap. We assume that each particle feels an equal 
mixture of Rashba and Dresselhaus spin-orbit coupling 
of strength k so and Raman coupling of strength Q. The 
system Hamiltonian H reads 

H = H sr I + ^ 1 Pi a y ) + } - C 1 ) 

i=i 

where 

N f.2 op. i 

H " = 'E~2^drf + 2^^ + ^ gS ( Xj ~ Xk) - (2) 

j =1 j<k 

In Eq. (|T]) , I is the 2x2 identity matrix, and a y ^ 
are spin-1/2 Pauli matrices of the jth particle, and pj 
is the momentum operator of the jth particle along the 
^-direction. Our goal in the following is to determine the 
eigenstates and eigenenergies E of H. 


To this end, we perform a unitary transformation and 
define H = WHU, where U = n^/ =1 exp(— ik so Xj(Ty^). 
The transformed Hamiltonian H reads 

H = H 0 i+^V R (3) 


with 


and 


Hq — Hsr 


NVkl 

2m 


N 

Vr = F exp (iksoXja^^a^ exp{-ik so Xja^). 
1=1 


( 4 ) 

( 5 ) 


The eigenvalues of H are the same as those of H while 
the eigenstates T of H are related to the eigenstates T of 
H by 'll = U Ik. Our strategy for solving the Sclirodinger 
equation ff’F = E^> is as follows: We first find the eigen¬ 
states of Hq and HqI and then account for the Raman 
coupling term Vr through either a matrix diagonalization 
or an effective low-energy Hamiltonian approach that is 
accurate to second order in fl. 

We denote the eigenstates and eigenenergies of H sr by 
4>n(x) and E n , where n collectively denotes the quan¬ 
tum numbers needed to label the states and x collec¬ 
tively denotes the N spatial degrees of freedom, x = 
(xi,X 2 ,- ■ -,Xn)- In general, the <fi n (x) are not known. 
Since H 0 is independent of spin, the eigenstates of H 0 I 
can be written as 


l/n,Si,S2 ,■■■,SN — 4 > n{x)\si, S2) ' ' '> (6) 

where the |si,S 2 ,- • -,sjv)j/ are eigenstates of the op¬ 
erator cjy = J2jLi a y\ i- e -> a< y\ s o)y = s i\ s i)y with 
Sj = ±1. The eigenenergies of these eigenstates are 
E n — NK 2 k1 0 /(2m). The eigenstates with the same <\> n (x ) 
but different spin configurations are degenerate. Thus 
one can form linear combinations of the '4 > n,si,s 2 ,---,s N such 
that the resulting eigenstates are eigenstates of the a 2 
operator, where r? 2 = &^) 2 - For non-vanishing 

H, If no longer commutes with <r y or a 2 , implying that 
the total spin and the total projection quantum numbers 
are no longer conserved. However, H still contains sym¬ 
metries. Specifically, H commutes with all Pjk (j, k = 
1 operators and the Y operator. Pjk exchanges 

both the spatial and spin coordinates of particles j and 
k. It follows that the eigenstates of H can be classified ac¬ 
cording to whether or not they change sign under the Pjk 
operator. Throughout this paper, we restrict ourselves to 
the N identical boson and N identical fermion sectors. 
The operator Y can be written as n jLiQj ® a x \ where 
Qj changes the sign of the spatial coordinate of parti¬ 
cle j. It follows that the eigenstates of H, which are 
also eigenstates of Y, always include an equal mixture of 
|Si, S 2 , • • •, S]sr)y and |si, S 2 , • • ■, s/v},,, where |sy) y denotes 
the state obtained by operating with onto \sj) y . Cor¬ 
respondingly, the expectation value of a y , and thus S y , 




3 


is always zero. We label the symmetries of the eigen¬ 
states by (a, b), where (a, b) = (1,1), (1, —1), (—1,1), and 
(-1,-1) indicate bosonic (a = 1) or fermionic symme¬ 
try (a = —1) with positive (b = 1) or negative (b = — 1) 
eigenvalue of Y. 

In what follows, we are interested in the expectation 
values of the local spins S x (x) and S z (x) in the x- and 
^-directions, 

N b 

S x{x) = ^2~a^ ) 8{x-Xj) ( 7 ) 

j =i 

and 

N 

s z (x) = 2 a z )s ( x - x i)- ( 8 ) 

j =i 

As already mentioned, we have ( S y (x )) = 0. To quantify 
the correlations of the system, we define the projector 
P\ M s |=mt 

P\M s \=m = E l Sl > s 2, ■ ' •, Sn) vv {si, S 2 , • ■ Sjv|, (9) 

| M s | =m 

where M s is the spin projection quantum number in the 

(/-direction, M s = S 1 +S 2 H-l-s/v- The expectation value 

°f P\m s \—m yields the probability to find the system in the 
space spanned by all the states with the same absolute 
value of M s . For example, for the three particle system, 
M s = ±1 and ±3. In this case, P\m s \=i and P\m s 1=3 
measure the configurations with M s = ±1 and M s = ±3, 
respectively. 

To diagonalize H directly, we start with the eigen¬ 
states of H 0 I [see Eq. ©]• By taking appropriate lin¬ 
ear combinations of il>n,ai,s 2 ,—,» jv , the Hamiltonian ma¬ 
trix for the different (a, b ) symmetry channels can be set 
up and diagonalized separately. The off-diagonal ma¬ 
trix elements of H contain spatial integrals of the form 


(f2/2) J2 x ,4 , n{x)4 , m{x)e ±2lkaoX:i dx, which is the Fourier 
transform with respect to the coordinate Xj of the prod¬ 
uct of two eigenstates of H 0 . As discussed in AnnendixlAl 
these integrals can be evaluated up to any precision nu¬ 
merically and in some cases can be calculated fully ana¬ 
lytically. 

Since we are primarily interested in the lowest eigenen- 
ergy and eigenstate of H for non-zero f2, we derive an ef¬ 
fective low-energy Hamiltonian, which describes much of 
the system dynamics accurately and is numerically more 
tractable. The idea is to divide the Hilbert s pac e spanned 
by H 0 I into two pieces called H L and H h [54], [55|. The 
low-energy space Hl contains the eigenstates with en¬ 
ergy less than a preset value and the high-energy space 
Hh contains the eigenstates with energy larger than this 
preset value. The eigenstates of the effective Hamilto¬ 
nian H e ff , which accounts for the Raman coupling term 
(Q/2) V/j perturbatively, are linear combinations of the 
unperturbed states of H 0 I that lie in Hl- We write 


where 

O’ _ o(0) , o(l) 1 o( 2 ) 

H eff - H eff + H eff + H eff’ 

(10) 


H eTf = E \^^,s 2 ,-,s N ) X 

Sl,S2,-",SjV 



S 1 ,S 2 ,---, S N 


(V’b,Sl,S2,- 

■■,SN\ A ljll\i’l ktS [ !S ' 2r .. !S 'J(^l kta ' 1 , 8 ' 2t ..., s ' N \ 

(11) 

with 

II 

(12) 


G |<M 
II 

^ -K 

(13) 


and 


Api* = ^ 


E 

// // / 
S 1 ) s 2 
hn £ Hh 


' v rW 


iXV’/, 


Eh — Eh„ 


,|Vk v R \\ n 

-— +-- 




,\Vr 


Ei k - E hn 


(14) 


( 2) 

An important point is that the second-order term A)./ 
runs over all states from the Hh Hilbert space. In 
essence, this leads to a renormalization of the matrix el¬ 
ements in the low-energy space Hl. 

To obtain an intuitive understanding of the effec¬ 
tive Hamiltonian approach, Figs, Oja) and QJb) show 
the single-particle momentum space potentials (p ± 
Hk so ) 2 / (2m) for k so aho = 2 and k so aho = 4, respectively. 


In this picture, the potentials centered at p = —Hk so 
and p = hk so are occupied by the spin-up and spin-down 
components of the wave function. The Raman coupling 
term introduces a mixing of the spin-up and spin-down 
potentials. Within the effective Hamiltonian approach, 
this is accounted for by H^jf and To estimate the 

relative importance of H^ff and H^f ■ we consider the 
g = 00 case (this case is discussed in detail in Sec. ns). 
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FIG. 1: Single-particle momentum space potentials (p±ftfc so ) 2 /(2m) (dashed lines) and non-interacting single-particle momen¬ 
tum space wave functions (solid lines) for (a) k so aho = 2 and (b) k so aho = 4, respectively. The potential and wavefunctions 
centered at p = —hk so are for the spin-up component and those centered at p = hk ao are for the spin-down component. For 
small k ao , the spin-up and spin-down single-particle wavefunctions in the Hilbert space Hl (here, states with E < 3.5 hui) 
have overlap, indicating that the first-order effective Hamiltonian dominates. For large k so , the spin-up and spin-down single 
particle wavefunctions in the Hilbert space Hl have essentially zero overlap, indicating that coupling to highly excited states 
is dominant, leading to a dominant second-order effective Hamiltonian. 


For g = oo, the 7V-particle wave functions are constructed 
from the non-interacting single-particle harmonic oscilla¬ 
tor states. For N = 4, e.g., the single-particle states 
shown in Figs. [H a ) and HP 3 ) contribute. The Hamil¬ 
tonian Hf.'f f couples states centered at p = —hk so and 
hk so . Visually, it is clear that the coupling (i.e., the 
overlap between the different single-particle states cen¬ 
tered at p = —hk so and Hk so ) decreases with increasing 
k so . As a consequence, the Hamiltonian term H x e ff may 

carry a higher “weight” than for sufficiently large 

( 2 ) 

k so despite the fact that is suppressed by the factor 
f l/(Hui) for H < Hco. More quantitatively, the TVth single¬ 
particle state is distributed in the momentum interval 
(±Hk so — hkF, ±hk so + /ffcp), where kpCLho ~ V2 N. Thus 
for k so > y/2N/aho , the coupling matrix elements enter¬ 
ing into are expected to be small while the coupling 

matrix elements entering into H x e ff are expected to have 
a comparatively large amplitude [the vertical arrow in 
Figs. [0(a) and[l](b) indicates the momentum region where 
the Nth single-particle state overlaps with a highly ex¬ 
cited single-particle state]. Our pictorial analysis is con¬ 
firmed by our quantitative calculations. Specifically, as 
detailed in Sec. lIVBl the spin structure changes drarnat- 
ically at a critical k so value, beyond which the term H y ^ 

dominates over the term • Analogous arguments can 


be made for finite g. 

Based on the effective Hamiltonian given in Eq. m, 
the next section calculates the spin structure for two- 
particle systems with arbitrary coupling constant g and 
compares the results with those obtained from the full 
(“brute force”) diagonalization. 


III. TWO-PARTICLE SYSTEM WITH 
ARBITRARY g 

For the two-particle system, the Schrodinger equation 
for H sr has been solved analytically for arbitrary two- 
body interaction strength in the literature [29]]. The 
eigenstates are cj> n (x ) = <h p (R)ip q (r) with eigenenergy 
E n = (p + 2q + where the center of mass co¬ 

ordinate R and the relative coordinate r are defined 
through R = (x\ + X2)/y/2 and r = (xi — a; 2 )/v / 2. 
The center of mass eigenstates < f’ p (i?) are the harmonic 
oscillator eigenstates with energy (p + 1/2 )hui, where 
p = 0,1, • • •. For states with even parity in the relative 
coordinate, the relative eigenstates p q {r) are given by 
N q U(-q 1 l/2,(r/a ho ) 2 )exp[-r 2 /(2al 0 )\, where a ho de¬ 
notes the harmonic oscillator length, 


aho = \ -, (15) 
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N q is a normalization constant and U is the confluent hy- culate the expectation values of S x (x) and S z (x), we find 
pergeometric function. In this even parity case, q denotes 
a non-integer quantum number, which is determined by 
the transcendental equation |29| 


(S b x (x)) = -C b x n b x (x) 


(18) 


2r(-g +1/2) 

r(-«) 


9 


V2huaho 


(16) 


For states with odd parity in the relative coordinate, 
the relative eigenstates <p q (r) are again harmonic oscil¬ 
lator eigenstates with eigenenergy (2 q + 1/2 )Hu>, where 
q takes half-integer values, i.e., q = 1/2, 3/2,- • •. Since 
the two-body d-function only acts at r = 0, the odd- 
parity eigenstates and eigenenergies are independent of 
g. The eigenenergies of the even-parity relative eigen¬ 
states, in contrast, change with g. Specifically, the en¬ 
ergy (2q+l/2)hui increases with increasing g. For infinite 
g , the ground state of the two-particle system is doubly 
degenerate, i.e., the relative energies of the even-parity 
state and the odd-parity state coincide, yielding a total 
energy of 2hui. This degeneracy plays, as we show be¬ 
low, an important role when the spin-orbit and Raman 
coupling strengths are non-zero. 

When k so and Q are non-zero, the center of mass and 
relative degrees of freedom are coupled. In order to de¬ 
scribe the interplay between the two-body interaction 
and the Raman and spin-orbit couplings, the states with 
spatial parts $o(R)ip qo (r) and $>o(R) l Pq 1 i r ) span the Hl 
space. Here tp qo and ip qi denote the energetically lowest- 
lying relative even-parity and odd-parity states, respec¬ 
tively. We choose the spatial basis functions to be real. 
Our motivation for choosing this Hl space is as follows. 
We want the resulting effective low-energy Hamiltonian 
to describe the ground state of the two-particle system 
with good accuracy for all g. For small g , e.g., the 
low-energy Hamiltonian constructed using &o(R)ip qo (r) 
would suffice since the energy of the state <S>o(R) ( Pq 1 (r) 
is close to huj higher in energy. Including this state in 
Hl changes the effective low-energy Hamiltonian and its 
resulting eigenenergies negligibly. For large <?, in con¬ 
trast, the state $o(R)<p qi (r) needs to be included in Hl 
since its energy is, as discussed above, nearly degenerate 
with the state <I>o (R)ip qo (r). Thus, the space Hl iden¬ 
tified above is the minimal space needed if an effective 
description of the ground state for all g (small and large) 
is sought. 

For two identical bosons, the ground state i/j b r of the 
effective Hamiltonian H e ff is 

^ b gr = ^MR)<P q o(r)(\n)y + \ll)y) 

+ ^<t> 0 (R)V qo (r)(\n)y + \Wy) 

+^o(R)9> q Ar)(\n)y-\lt)y). (17) 

The coefficients C\-C^ are obtained by diagonalizing the 
effective low-energy Hamiltonian H e ff. Using i/j b r to cal- 


and 


where 


(S b (x)) = g C> b z (x), 


c b = c b (c b r + ( C b yc b , 


(19) 


( 20 ) 


/ oo 

$l(R)<p 2 qo (r)[6(x-x 1 )+6(x-x 2 )]dx, (21) 

-oo 


C b z = i[C b {C b r - (C b )*C b ], 


( 22 ) 


and 


/ OO 

$o (-ft V® (O^gi ( r ) [5{x-xi)-5{x- x 2 )\dx. 

-OO 

(23) 

For two identical fermions, the ground state of the ef¬ 
fective Hamiltonian is 

nf 

^ gr = -±MR)‘P q Ar)(\n)y + \tt)y) 

nf 

+ ^ Q (R) lpqi (r)(\n)y + \lt)y) 

+ ^ 0 (R) ipqo ( r )(\n)y-\Wy). 


(24) 


It can be seen that the structure of ipf gr is very similar 


to that of ip b gr , Eq. m The only differences are that 
the superscript b is replaced by /, that ip qo is replaced 
by <p qi (twice) and that ip qi is replaced by ip qo (once). It 
follows that the expressions that describe the spin struc¬ 
ture of the fermionic system are nearly identical to those 
that describe the spin structure of the bosonic system. 
Specifically, Eqs. (HI-® apply to the fermionic system 
if the superscript b is replaced by /, ip qo is replaced by 
ip qi , and ip qi is replaced by ip qo . 

We refer to Cll^ and C b J'^ as the spin structure coef¬ 
ficients and to njf{x) and nj^\x) as the spin structure 
densities. Note that these densities can be negative; a 
negative density indicates that the spin points in the neg¬ 
ative direction. Importantly, the spin structure densities 
depend on g but are independent of k so and Q. The spin 
structure coefficients, in contrast, depend on k SOl fl and 
g. To obtain the spin structure within the effective low- 
energy Hamiltonian approach, the spin structure density 
(see Fig. [3]) gets multiplied by the spin structure coeffi¬ 
cients, shown in Fig. [2] as a function of k so for various g 
and fixed f l (fl = Huj/2). 






0 1 2 3 4 5 

kso^ho 


FIG. 2: Spin structure coefficients (a) C x , (b) C b , (c) C 1 /, 
and (d) C( for the N = 2 ground state obtained within 
the effective Hamiltonian approach for 12 = ftuj/2. The 
spin structure coefficients in (a) and (b) are for the sys¬ 
tem consisting of two identical bosons. The spin struc¬ 
ture coefficients in (c) and (d) are for the system consist¬ 
ing of two identical fermions. Dashed, dot-dashed, dot-dot- 
dashed, and solid lines show the spin structure coefficients 
for g = y/2hwaho , 3 \/2hu>aho, 1 5V2hu)aho and oo, respectively. 
For infinitely large g (solid lines), the spin structure coeffi¬ 
cients display a “sharp spike” at k so a,ho ~ 3.4; this is the 
result of an avoided crossing between the ground and first 



FIG. 3: Spin structure densities for the N = 2 ground state 
for 12 = frw/2 and varying g. Panels (a) and (b) show the spin 
structure densities n b (x ) and n b (x), respectively. Dashed, 
dot-dashed, dot-dot-dashed, and solid lines show the spin 
densities for g = y/2hcuaho, 3y/2hujaho, 15y/2hujaho, and oo, 
respectively. The spin structure density is independent of 
k so - The spin structure density nl(x) is independent of g and 
shown by the circles in (a). The spin structure density n{{x) 
coincides with n b (x) for all g. 


In general, the spin structures of the bosonic and 
fermionic systems differ. For a weak two-body inter¬ 
action (small positive g) and Raman coupling strength 
f 1 = hoj/2 with vanishing spin-orbit coupling strength, 
the coupling between the singlet and the triplet states 
vanishes. Since the relative even parity state has a 
lower energy than the lowest relative odd-parity state 
for g < oo, the ground states for two identical bosons 
and fermions are triplet and singlet states, respectively. 
Correspondingly, the coefficient in Eq. 03 and the 



FIG. 4: Benchmarking the effective Hamiltonian approach for 
the N = 2 ground state for Q = hu/2 and g = oo. Panels (a) 
and (b) show (S x (x)} and (S z (x)) for k so aho = 1/5; panels 
(c) and (d) show (S x (x)} and (S z (x)} for ksoaho = 4. The 
solid lines show the spin structure obtained from the exact 
diagonalization while the dashed lines show the spin structure 
obtained within the effective Hamiltonian approach. On the 
scale shown, the solid and dashed lines nearly coincide. 


coefficients C( and C( in Eq. (EH) vanish, yielding 
(^0)> = {Sf(x)) = (S((x)) = 0 and (S b {x)) ^ 0. 
For a non-vanishing spin-orbit coupling strength, the sin¬ 
glet and triplet states are coupled, leading to non-zero 
(S b (x)), (S£(x)), and (S£(x)). With increasing g , the 
energy difference between the singlet and triplet states 
for k so = 0 decreases. As a consequence, for a fixed and 
relatively small k so , the spin structures for two identical 
bosons and two identical fermions are more similar for 
larger g than for smaller g [see Figs. El[a) and Etc) for 
k so (iho ;$ 1.5]. However, the coupling between the sin¬ 
glet and triplet states is weakened with increasing k so . 
Specifically, for fixed g 1 the spin structures for bosons 
and fermions differ more as k so increases [see Figs. [2ja) 
and Etc) for k so aho > 2]. For a larger k so , a larger g 
is needed to get the triplet (singlet) state mixed signif¬ 
icantly into the ground state for two identical fermions 
(bosons), resulting in similar spin structures for bosons 
and fermions. For an infinitely large g 1 the relative even 
and odd parity states have the same energy and the spin 
structures for bosonic and fermionic systems are identi¬ 
cal for all k so and 12 (see the solid lines in Fig. EJ. This 
effect can be attributed to the Bose-Fermi duality (see 
next section for more details). 

Our discussion so far has been based on the effective 
Hamiltonian approach. To benchmark this approach, we 
compare the effective Hamiltonian results with those ob¬ 
tained from the exact diagonalization for various g , k so , 
and 1 1 combinations. The agreement is good for all cases 
considered. As an example, Fig. El compares the two- 
particle spin structure for infinite g and 12 = hio/ 2. Fig¬ 
ures E[a) and [D(b) are for small k so ( k so aho <C 1) and 
Figs. Etc) and[D(d) are for large k so (k so aho 1)- The 
agreement between the effective Hamiltonian approach 
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results (solid lines) and the exact diagonalization ap¬ 
proach results (dashed lines) is very convincing. 


X 2 < ■ ■ ■ < Xn- Within this subspace, the evaluation of 
the Hamiltonian matrix elements involve integrals of the 
form 


IV. N-PARTICLE SYSTEM WITH g = oo 
A. Formulation 


V 


ni ,n 2 .. .njv 


D ( n l 


5 '*2» 


, UN 


)X 


772 , * * *, TIn^)®x\<x-2---<xn^ S ° ^dx. 


(26) 


When the two-body coupling constant g is infinitely 
large, the particles cannot pass through each other. 
In the absence of the spin-orbit and Raman coupling 
terms, the atomic gas behaves like a Tonks-Girardeau 
gas [30, [63] ■ The Tonks-Girardeau gas has a large de¬ 
generacy and bosons fermionize {Hi, j63|, 0|. The fact 
that the particles cannot pass through each other im¬ 
plies that the particles can be ordered. Since there are 
N\ ways to order the particles, the degeneracy of each 
eigenstate of H sr with g = oo is N\ if the particle ex¬ 
change symmetry is not being enforced. We thus pursue 
an approach where we first determine the eigenstates of 
H sr I for a fixed particle ordering. The resulting eigen¬ 
states are then used either to calculate the eigenstates 
and eigenenergies of H through exact diagonalization or 
to calculate the eigenstates and eigenenergies of H e f / by 
identifying an appropriate subspace Hl. The fact that 
the particles can be ordered also allows us to derive an 
effective spin Hamiltonian that is independent of 

the spatial coordinates. 

For infinite p, the eigenstates of H sr can be written 
as [3lj ] 


<f>n 1 ,n 2 ,...,n N (x) — 77,2, ' ' njv)0xj 1 <Xj 2 <■■■< Xj N , 

(25) 

where D{n\,ri 2 , • • -,njv) denotes the Slater determi¬ 
nant constructed from N one-dimensional harmonic os¬ 
cillator eigenfunctions with quantum numbers ni,772,- • 
■,riN (n i ^ 772 7^ • • • 7^ 77jv)- The sector function 
0x jl<Xj2 <-<x jN is 1 for x h < Xj 2 < ■ ■ ■ < x jN and 
zero otherwise. For the ground state, e.g., we have 
77i = 0 ,772 = 1, • • -jIT-n = N — 1. To construct eigen¬ 
states of H sr I , we include—as before—the spin part 
|si,S 2 j - • ■, S]y)y. It is important to realize that there is 
no coupling between the eigenstates for different particle 
orderings, i.e., the Raman coupling term (: Q/2)Vr does 
not couple states with different particle orderings. This 
implies that the Hilbert space can be divided into N\ in¬ 
dependent subspaces. A state in a given subspace can be 
mapped onto a state in a different subspace with the same 
eigenenergy through the application of one or more per¬ 
mutation operators. For example, for the N = 3 system, 
the particle ordering X\ < X 2 < £3 can be changed to the 
ordering X 2 < X\ < X 3 through the application of the P \2 
operator. This property is used, after solving the problem 
in one of the N\ subspaces, to construct fully symmet¬ 
ric bosonic or fully anti-symmetric fermionic eigenstates 
from the non-symmetrized eigenstates that span one of 
the N\ distinct Hilbert spaces. 

Without loss of generality, we discuss the ordering X\ < 


The evaluation of this integral for any N is de¬ 
tailed in Appendix [A] Once the Hamiltonian ma¬ 
trix has been diagonalized, we determine the fully 
symmetric/anti-symmetric eigenstates by applying the 
IV-particle symmetrizer/anti-symmetrizer. 

Since we consider a fixed particle ordering, we define 
the sector Hamiltonian H Xji<Xj2< ... <x . N , 


Hx jl <x j2 < - <xj N — H0x jl< x j2< - <XjN , ( 27 ) 

which is non-zero only for the particle ordering 
Xj 1 < Xj 2 < • • • < Xj N . From the defi¬ 

nitions of H Xji<Xj2< ... <XjN and Y, it follows that 
YH X ji<x j2 < - < x jN = i i- e -> the 1 

operator does not commute with the sector Hamil¬ 
tonian. We can, however, define a sector operator 
y Xjl<Xj2< - <XjN for the particle ordering x jl < x n < 
• • • < Xj N , which commutes with the sector Hamilto¬ 
nian with the same particle ordering. The eigenstates 
of the Y operator are then obtained by symmetrizing or 
anti-symmetrizing the eigenstates of the Y Xj ^ < Xj2 <---<x 3N 
operator. The key idea is that the particle ordering after 
a change of the spatial and spin coordinates can be “re¬ 
stored” by exchanging the coordinates of the particles. 
Specifically, Y Xji<Xj2< ... <XjN is defined through 


Y, 


Xj, <Xn 0 <-"<Xn 


= (±i)M 


P p . p y 

^3l3N r 323N-l ^3 [ N ] 0 NJrl _ [ N ] 1 J 


(28) 


where [iV/2] denotes the integer part of N/ 2, and the 
plus and minus signs apply to identical bosons and 
fermions, respectively. It can be proven readily that 
Y Xji<Xj2 <...< XjN commutes with , 1A . To 

show that the eigenstates of the Y operator are obtained 
by symmetrizing or anti-symmetrizing the eigenstates of 
Y Xjl<Xj3< ... <XjN , let <x !2 < — x :l> . denote a wave func¬ 

tion in the subspace with the ordering < Xj 2 < ■ ■ ■ < 
Xj N with the property 


Y 


<Xn 0 <‘<Xj 


r V’ayi 

= 


<Xn 0 < - <Xn 


< x io 


(29) 


Taking advantage of 


S PnmPj 2jN -i '' ' Pj,K]0 N+ i-[K] ~ S i (30) 


S a {~lp ] Pno N Pj 2 m-^ -P> 


■INI N = S a , (31) 


and the fact that Y commutes with S s ^ , where S 
denotes the JV-particle symmetrizer (anti-symmetrizer), 
we obtain 


S s(a) Y 


Xji <Xj 0 <---<Xj 


YS s ^ip Xji 


< ^ x io <••• < C.Xj 


< Xjn <‘”< Xj 


(32) 


. <Xj 0 <---<Xi 


Eq. (13211 leads to 
= bS s( ' a ' > ip Xji<Xj2< ... <XjN , 
is an eigenstate of Y with 


Using Eq. ([29]) . 

1 S s ^i/j Xji<Xj2< ... <x 

i.e., S< a ^ Xji 
eigenvalue b. 

To construct the Hamiltonian H e ff 1 we proceed simi¬ 
larly, i.e., we also work with a particular particle ordering. 
The space Hl is spanned by the states D(0,1, • • -,N — 
1 )©xi<x 2 <---<xjvI s 1 ) s 2 , ■ • •, s N ) y , where the spin function 
can take 2 N different arrangements. Correspondingly, 
the space Hh is spanned by all other unperturbed eigen¬ 
states with the ordering X\ < X 2 < ••• < Xn ■ The Hamil¬ 
tonian matrix for H e ff is constructed and diagonalized, 
and states with good symmetry are obtained following 
the same steps as discussed above. 

We now discuss the construction of an effective spin 
Hamiltonian. Since the particles can be ordered in N\ 
distinct ways, the ground state is iV!-fold degenerate. 
We integrate out the spatial degrees of freedom of the 
effective Hamiltonian H e ff , yielding a spin Hamiltonian 
H s e P f l p that depends only on the spin degrees of freedom. 
Specifically, since the Hilbert space Hl contains exactly 
one spatial wavefunction with ordering X\ < X 2 < • • • < 
Xn, we define H s e pl p through H*jJ l = f^° OC) D*( 0,1,- • 
•,7V- l)H e f f D(0, 1, - - TV — l)0xi <x 2 <-<*i ydx. Using 
that the spins are also ordered, H s e pl p can be compactly 
written as 


KTF = [Eo- 


NVkl 

2 m 


n 


N 


EE 


3 = 1 


~ 2 hui E a J M i k<Tk + Ah,,, EE + 


j<k 


where 


and 


= (< 


4 hio 


r0‘) n-(l) 


Bj = (B«,BW). 


(33) 


(34) 


(35) 


Eq is the ground state energy of H sr , E 0 = HuiN 2 / 2 , and 
M 3 i- is a 2 x 2 matrix (see below). The first, second, and 
third terms on the right hand side of Eq. (1331) come from 
and respectively. The fourth term on 

the right hand side of Eq. (1551) also comes from 7?e//> 
it accounts for the (spin-independent) onsite interaction, 
with ajj and bjj defined below in Eqs. m and m - We 
find 


= Re{Vi 


0,1,2, —,JV-l) 
0,l,2,---,iV—1 


(36) 


and 


0,1,2, —,iv-1 

The matrix Mjk can be written as 


Mjk = 


Qj k Cjk 
Ckj bj k 


( 37 ) 


(38) 


where 


Ujk 

hui 


bjk_ 

Hlo 


E 

(m,n 2 ,n3,-,niv) 
#(0,1,2,-,JV-1) 


-^ e (E,l,2, —,iV-l)^ e (E,l,2,—,AT-l) 
Til ,ri2 ,--,TlN Tli,ri2Y" i n N 


Ef) — Er 


ni,n 2 ,’”,TiN 


E 

(n 1 ,n 2 ,n 3 ,---,n N ) 
/(0,1,2,-,JV—1) 


(39) 

^(®S,l,2,..,lV_l)^(®5 ) l,2,.. 1 lV_l) 

ni,n 2 Y‘‘ i n N Tii,n 2 ,---,riN 


Eq — E, 


ni,n 2 ,---,nN 


C jk 

hui 


and 


c kj 

hui 


E 


(40) 

^ e (^0,l,2,—,JV-l)-^ m (®0,l,2, —,JV-l) 
ni ,ri 2 7^1,7i2,---,njv 


(ni,n 2 ,ri3,-",niv) 

^(0,1,2,-,JV-1) 


En — E n 


(41) 


E 

(ni,n 2 ,n 3 ,---,n N ) 

^(0,1,2,-,1V-1) 


I' rn ('Db t l, 2 ,—,N-l)R e (E > 0 ,l, 2 ,—,N-l) 

ni ,n 2 ,---,TlN Til,7l 2 ,--,TIN 


Eo — E, 


n i ,n 2 ,---,riN 


(42) 

The second and third terms on the right hand side of 
Eq. (1551) correspond to a single spin term (this is the first- 
order term) and a spin-spin term, respectively. For small 
k so , the first-order term dominates. The spin at each slot 
j follows the effective magnetic B, field and the system 
has a rotational spin structure [2l|. In this regime, the 
spin correlations are very weak. The eigenstates of 
include equal weights of each spin state, namely, P\m s \ is 
proportional to the number of spin states that have the 
same absolute value of M s (see Table [T|. For example, 
for N = 3, two states have \M S \ = 0 and six states have 
\M S \ = 1 , which yields that P\m s \=o = 1/4 and P\m b \=i = 
3/4. 

For fairly large k so , the second-order term dominates 
over the first-order term, i.e., the spin correlations are 
very strong and the spin structure is non-trivial (see also 
discussion around Fig. [TJ. For large k so , we find nu¬ 
merically that the nearest neighbor spin-spin interactions 
dominate. For example, |a.i 2 1 is much larger than |< 3 , 13 1. 
It should be noted, though, that the summands enter¬ 
ing into ai 2 and 013 are of roughly the same order of 
magnitude. Moreover, we find that the nearest-neighbor 
coefficients a^jk),bijk),cijk) and Ctkj) are, for fixed k so , 
approximately equal; here, the subscript (jk) indicates 
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that j and k are related via j = k — 1. We have checked 
this for IV < 4. For large N, we cannot rule out that 
the values of the coefficients depend on the slot, possibly 
yielding more complicated spin structures than discussed 
in this work. For a^jk) = = C(jfc) = the matrix 

Mjk simplifies dramatically, 


M (jk) 



-1 

1 


(43) 


for j = k — 1 and Mjk ~ 0 for j ^ k — 1. In Eq. (l43l) . 
A is a negative dimensionless constant. In this case, the 
last term on the right hand side of Eq. (1551) reduces to 


n 2 

2 hbj 


Y. drjMjkOk 

j<k 



A y [°j • °k - y x j) ^i k) - 

(jk) 


(44) 


The first term in square brackets on the right hand side of 
Eq. (l44l) corresponds to the usual Heisenberg exchange 
term. The term in round brackets on the right hand 
side of Eq. (1511) is readily indentified as the y-component 
of the cross product between two three-dimensional spin 
vectors. This term is the one-dimensional analog of the 
anisotropic Dzyaloshinskii-Moriya exchange term [60l — 

m. 

To obtain the spin structure of this approximate effec¬ 
tive spin Hamiltonian, we rewrite Eq. (1441) as 


H 2 

2 Huj 


yy MjkOk 

j<k 


n 2 

2 tux 


Ayi 2(1 - i)o ( l ] o ik) + 2(1 + ijcx^V* 0 ], 
( jk ) 


(45) 


where <j±^ = (cxi^ p icr%)/2. If we neglect the first- 
order term and use the right hand side of Eq. (1551) in 
Eq. (1551) . we find that commutes with o y . More¬ 
over, this approximate Hamiltonian also has time 

reversal symmetry, i.e., it commutes with io y K, where 
K changes the quantity that it acts on into the com¬ 
plex conjugate of that quantity. However, io y K and o y 
do not commute. Using additionally the property that 
K\s) y = |s) y , one can show that the eigenstates with M s 
and — M s of the approximate effective spin Hamiltonian 
have the same energy provided \M S \ > 0, i.e., the states 
with | M s | > 0 are two-fold degenerate. The first-order 
term breaks the degeneracy of the eigenstates with M s 
and —M s . As a result, the eigenstates of the full Hamil¬ 
tonian -Hfy™ (including zero-, first-, and second-order 
terms) are, for k so —> 00 , approximately superpositions 
of states that have the same absolute value of the M s 
quantum number. The 0 ^ 0 ^ and <7^%% terms cor¬ 
respond to nearest neighbor spin hopping. These terms 
lead to a lowering of the energy. The more possibility 
for the nearest neighbor spin hopping a state has, the 


N 

kso — y 0 

kso — t 00 

2 

P 0 = 1/2, P 2 = 1/2 

Po = 1, P 2 = 0 

3 

Pi = 3/4, P 3 = 1/4 

Pi = 1, P 3 = 0 

4 

Po = 3/8, P 2 = 1 / 2 , P 4 = 1/8 

Po = 1, P 2 = 0, P 4 = 0 


TABLE I: Spin correlations for the ground state in the limit¬ 
ing cases kso —» 0 and k BO —» 00 for various N. The probability 
P\m s \ that the ground state has the absolute value of M s is 
reported. For small k so , all spin states are equally weighted, 
which means that P\m b \ is proportional to the number of spin 
states that have the same \M S \. For large k 30 , the spin states 
are not equally weighted. In this case, the ground state con¬ 
tains only spin states with the minimum allowed \M a \. 


lower the energy associated with that state is. For ex¬ 
ample, the states | f4% and | 44% are “connected” via 
nearest neighbor hoppings while the states | f1% an d 
| 44% are not connected with each other or with | 14% 
or | 44% via nearest neighbor hopping. Correspondingly, 
the N = 2 ground state is a linear combination of the 
| t4% and | 44% states. For N = 3, e.g., the state | tt4% 
is connected to | 144% via ° 2 _) ° 3 + ' 1 and to | 444% via 
txJ _ % 3 + ' ) while the state | 444% is connected to | 444% 

via o' 2 + % 3 ~' 1 and to | 144% via cr^o^ K Correspond¬ 
ingly, the N = 3 ground state is a linear combination of 
all \M S \ = 1 states. Table Q] shows the values of P\M s \=m 
in the limits k so —> 0 and k so —> 00 for N = 2 — 4. 

Figure O shows the probability to find the system in a 
state with a given \M S \ as a function of k so for the two- 
, three-, and four-particle systems with infinitely large 
g (see the next subsection for the calculational details). 
For large k so , P|m 3 |=o = 1 f° r the two-particle system, 
P\m s \=i = 1 for the three-particle system, and P\m s \=o = 
1 for the four-particle system. 


B. Application to systems with N = 2 — 4 

This section evaluates the expectation values of the 
spin operators defined in Eqs. 0-© for the three- and 
four-particle systems with infinite g. For N = 2, the 
expectation values of S x (x) and S z (x) for infinitely large 
g have been calculated in Sec. HTT1 For N = 3, the non- 
symmetrized ground state wave function of the effective 
Hamiltonian H e ff is 


^(0,1,2% 

Ygr — ™X 1 <X2<X3 x 

Cl (I 144% - I iU)y) + C 2 (I m>„ - I 144%) 


+%(l 144% - 1 444%) + c 4 (1 444% - I W1%) • (46) 


The coefficients C 1 -C 4 are obtained by diagonalizing the 
effective low-energy Hamiltonian H e ff. Using this wave 
function, the expectation values of S x (x) and S z (x) for 














FIG. 5: Expectation value of P\M s \=m for (a) TV = 2, (b) TV = 
3, and (c) TV = 4 with g = oo and Q = fko/2 as a function of 
k 30 ■ (a) The solid and dashed lines show P\m b \=o and P\m„\= 2 , 
respectively, (b) The solid and dashed lines show P\m b \=i and 
P | m b |= 3 ) respectively, (c) The solid, dashed, and dot-dashed 
lines show P\m b \=o, P\m b \= 2 and P|m s |= 4 , respectively. 


the three-particle system are 

(S x (x)) = ^[C lx n lx (x) + C 2x n 2x (x)\ (47) 

and 

(Sz(x)) = ^ C z n z (x ). (48) 

The coefficients C\ x ,C 2x , and C z , which depend on the 
coefficients C\ — C 4 in Eq. (l46l) . and the expressions for 
the density functions n\ x {x), n 2x [x), and n z (x) are given 
in Appendix [B] Figure [6] shows the spin coefficients and 
spin densities for TV = 3 with infinite g for different k so . 
For fixed x, ( S x (x )) and (S z (x)} oscillate with k so for 
k so < kgg , where a/jofcgo « 3.7. The oscillations disappear 
for k so greater than 

For TV = 4, the non-symmetrized ground state wave 
function of the effective Hamiltonian H e f f is 


i/j gr =D{ 0,1,2,3)- 


v^L 


Cl (| tttt) y + I UU)y) 


+c 2 (i tm), +1 wu) s ) + c 3 (i mt)y ■ 
+C 4 (I nn)y + I UU)y) + C 5 (I 


Uit )y) 


+c 8 (i mt)y +1 mi)v) + Cei tm)y + c 7 \ nn), 


+C9\un)v+c 10 \mf) v > ©a:i <x 2 <x 3 <.) 


The coefficients Ci-Cio are obtained by diagonalizing the 
effective low-energy Hamiltonian H e ff- Using this wave 



FIG. 6: Spin structure for the three-particle system with g = 
00 and fl = fojj/2 as a function of k so ■ (a) The solid and 
dashed lines show the spin structure coefficients C\ x and C 2x , 
respectively, (b) The solid and dashed lines show the spin 
structure coefficient C z . (c) The solid and dashed lines show 
the spin structure densities n\ x {x) and n 2x (x), respectively, 
(d) The solid line shows the spin structure density n z (x). 


function, the expectation values of S x (x) and S z (x) for 
the four-particle system are 

(S x (x)) = ^\Ci x iii x {x) + C 2x n 2x (x)] (50) 

and 

(S z (x)) = ^\C lz n lz {x) + C 2z n 2z (x)]. (51) 

Expressions for the coefficients Ci x ,C 2x ,C\ z , and C 2z 
and the density functions ni x (x),n 2x (x),ni z (x), and 
n 2z {x) are given in Appendix iBl Figures [71(a) and [71(b) 
show that the spin structure coefficients go to approxi¬ 
mately zero at ah o k^ r 0 « 4.3. 

To visualize the spin structure, Figs.ISlfHllshow the spin 
vector ((S x {x)), (S z (x))) at each “slot” as a function of 
k so for N = 2 — 4 with g = 00 . In each figure, the top 
to bottom panels correspond to the leftmost to the right¬ 
most slot. For k so ;$ the spin vector rotates with 
increasing k so ■ For k so > the spin-spin interaction 
term dominates. For N = 2 and TV = 4, the magnitude 
of the spin vector at each slot is approximately zero for 
k so ^ fcgo- For TV = 3, in contrast, the magnitude of the 
spin vector at each slot is fixed and the orientation of 
the spin vector is, to a good approximation, independent 
of k so . Taking the ground state of the approximate ef¬ 
fective Hamiltonian to calculate (S x (x)) and (S z (x)), we 
find ( S x (x)) ~ (, S z (x )) ~ —h/2 for the leftmost slot and 
(, S x (x )) ~ —(S z (x)) ~ — h/2 for the rightmost slot for 
k so ^ fcgo- This behavior is a signature of the spin-spin 
correlations. For even TV systems, the ground state in 
the large k so limit is a superposition of spin states with 
\M S \ = 0. Two states |si, s 2 , ■■■, sn) v and |s lf s 2 , s N ) y , 
both with M s = 0, contain an even number of Sj and Sj 






















FIG. 7: Spin structure for the four-particle system with 
g = oo and ft = hx/2 as a function of k so . (a) The solid 
and dashed lines show the spin structure coefficients C\ x and 
C 2 x■ (b) The solid and dashed lines show the spin structure 
coefficients C\ z and C 2z ■ (c) The solid and dashed lines show 
the spin structure densities n\ x {x) and ri 2 x(x). (d) The solid 
and dashed lines show the spin structure densities ni z (x) and 
n 2z (x). 



k £L 

so no 



FIG. 9: The spin structure for three identical particles with 
g = o o and fl = hco/2. Panel (a) shows the spin vector at the 
leftmost slot, panel (b) shows the spin vector at the middle 
slot, and panel (c) shows the spin vector at the rightmost 
slot. For k ao < kgo « 3.7/a ho, the spin vector at each slot 
follows the local effective B field. The spin vector rotates 
with increasing k ao . For k ao > k 3 „, the spin-spin interaction 
is dominant, resulting in a spin vector with approximately 
constant magnitude and orientation at each slot. The inset 
in (b) replots the spin vector using a different orientation of 
the coordinate system. 


((S x (x)), ( S z (x )}) vanishes for even N systems with large 
k so while that for odd N systems is finite and approxi¬ 
mately constant. 


V. CONCLUSION 


FIG. 8: The spin structure for two identical particles with 
g = oo and fl = hw/2. Panel (a) shows the spin vector 
at the leftmost slot while panel (b) shows the spin vector 
at the rightmost slot. For k so < k 3 „ ~ ‘iA/aho , the spin 
vector at each slot follows the local effective B held. The spin 
vector rotates with increasing k so ■ For k so > k 3 „, the spin- 
spin interaction is dominant, resulting in an approximately 
vanishing spin vector at each slot. 

for which Sj ^ s . For N = 2, e.g., to get the state | f4% 
from the state ] 44%) two spin flips are needed. For 
odd N systems, the ground state in the large k so limit 
is a superposition of spin states with \M S \ = 1. Two 
states |si,S 2 )- • • ,SN)y and |s 1 ,s 2 ,- • -,s N ) y , both with 
| M s | = 1 , contain an odd number of Sj and Sj for which 
Sj ^ Sj. For N = 3, e.g., to get the state | ft4% from the 
state | t44%, one spin flip is needed. The non-vanishing 
(S x {x)) and ( S z (x )) arise from a superposition of states 
that contain N — 1 identical spins (Sj = Sj ) and one pair 
of opposite spins (Sj ^ Sj). As a result, the spin vector 


Spin Hamiltonian play an important role in under¬ 
standing material properties such as the transition from 
ferromagnetic to anti-ferromagnetic order. Much effort 
has gone into identifying clean model systems with which 
to emulate spin dynamics. Notable examples include 
trapped ion systems, ranging from small two- or three- 
ion chains [bdl . to large two-dimensional ion crys¬ 
tals m HI , and ultracold atoms [H, [T?], [5614591 [gqI [7Cj ■ 
This paper considered ultracold harmonically trapped 
one-dimensional atoms with infinitely large two-body 
contact interactions subject to spin-orbit and Raman 
couplings and derived an effective low-energy spin Hamil¬ 
tonian that is accurate to second order in the Raman cou¬ 
pling strength f2. It was shown explicitly for TV = 2 — 4 
particles that tuning the spin-orbit coupling strength 
from a regime where the spins independently follow the 
effective external magnetic field generated by the Raman 
coupling to a regime where the spin dynamics is governed 
by the nearest neighbor spin-spin interactions is feasible. 
While the examples presented are for small particle num¬ 
bers, Appendix O derived a number of identities appli¬ 
cable to systems of arbitrary size that should facilitate 
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FIG. 10: The spin structure for four identical particles with 
g = oo and fl = Tujj/2. Panel (a) shows the spin vector at 
the leftmost slot. Panel (b) shows the spin vector at the 
second slot from the left. Panel (c) shows the spin vector at 
the third slot from the left. Panel (d) shows the spin vector 
at the rightmost slot. For k so < ~ 4.3 /aho, the spin 

vector at each slot follows the local effective B field. The spin 
vector rotates with increasing k so . For k so > the spin- 
spin interaction is dominant, resulting in an approximately 
vanishing spin vector at each slot. 


extensions to larger N and that may find applications 
in calculating the momentum distribution or related ob¬ 
servables for strongly-interacting one-dinrensional atomic 
gases without spin-orbit coupling. 

The relevance of the effective spin Hamiltonian derived 
in this paper is two-fold, (i) As already alluded to above, 
this paper introduced a new means to realize a tunable 
effective spin Hamiltonian, (ii) The spin Hamiltonian 
language provides a physically transparent means to un¬ 
derstand the intricate dynamics of strongly correlated 
spin-orbit coupled ultracold atomic systems. 

Throughout this paper, explicit calculations were per¬ 
formed for the Raman coupling strength fl = hoj/2. For 
this coupling strength, the energy difference for k so > k^ r 0 
between the ground state and the first excited state is of 
the order of 10~ 4 huj for N = 2 — 4 particles. This small 
energy difference may make it challenging to experimen¬ 
tally occupy the ground state. Since the second-order 

I 


term is proportional to SI 2 /(Hui), the critical k^ r 0 decreases 
with increasing fl. For Q = 4hui and N = 3, e.g., we find 
kso a ho ~ 3.5 and an energy splitting between the ground 
state and the first excited state with the same symmetry 
of the order of 10 ~ 2 hoj. This energy splitting appears 
much more tractable experimentally. 

An important question is then whether the effective 
second-order low-energy spin Hamiltonian is applicable 
for such large H. To investigate this question, we com¬ 
pared the results obtained by solving the second-order 
effective Hamiltonian with the results obtained from the 
diagonalization of the full Hamiltonian. The agreement 
is very good for k so > (we excluded the regime 
k so suggesting that the third-order term, which 

is proportional to H 3 /(feu) 2 , is suppressed. Naively, one 
expects the third-order term to be enhanced by kl/(hui) 
compared to the second-order term for H > fno , making 
a perturbative approach meaningless. However, the full 
perturbative expression also contains the matrix elements 
and energy denominator. It turns out that the product 
of three matrix elements is highly suppressed compared 
to the product of two matrix elements. We conclude that 
the second-order effective spin Hamiltonian capture the 
physics, including the spin structure up to, at first sight, 
surprisingly large fl/(fiu) for relatively large k so . 

A key ingredient that went into deriving the effective 
low-energy spin Hamiltonian is that the particles in one¬ 
dimensional space can be ordered if the two-body cou¬ 
pling constant g is infinitely large. For finite g , this is 
not the case, i.e., particles are allowed to pass through 
each other. In this case, a low-energy Hamiltonian that 
depends on the spatial and spin degrees of freedom was 
derived (in fact, the effective spin Hamiltonian for g = oo 
was derived by taking this Hamiltonian and integrat¬ 
ing out the spatial degrees of freedom). The low-energy 
Hamiltonian was tested for two particles and shown to re¬ 
produce the full Hamiltonian dynamics well. Moreover, 
it was shown to provide a powerful theoretical frame¬ 
work within which to interpret the full Hamiltonian re¬ 
sults. We believe that the formalism can be applied to 
larger one-dimensional system and extended to higher¬ 
dimensional systems. 
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Appendix A: Calculation of involved integrals 

This section contains the evaluation of the integral given in Eq. (25D. We denote the nth harmonic oscillator 
eigenstate by ip n {x), <Pn(x) = N n H n (x/ah 0 )e~ x where N n = 1/^/nn\2 n aho is the normalization constant 
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and H n {x ) the nth Hermite polynomial. Throughout this appendix, we set aho = 1, i-e., we work with dimensionless 
spatial coordinates. Expanding the Slater determinant D(ni,n 2 , • • •, n jv), we have 

D(ni,n 2 ,-■ -,n N ) = (~ 1 ) Ppi ' P 2 '"' PNll iLi l fn Pl { X i), (Al) 

Pl,P2,-",PN 

where pi,p 2 ,- --,Pn denotes a permutation of 1, 2, • • •, TV and Pp 1 ,p 2 ,---,PN the number of permutations needed to obtain 
the order pi,p 2l - ■ -,pn from the ordinary order 1,2,- • -,N. The sum in Eq. EH) contains N\ terms. Note that 
since the eigenstates </> raiin 2 ) ... jnw (x) in Eq. E5l) are only non-zero for a particular particle ordering, we do not need 
a prefactor of l/\/iVT in front of the Slater determinant to normalize the eigenstates. Equation (l26l) contains two 
different Slater determinants (D functions), one with arguments n\, ■ ■ ■, n/v and the other with arguments n l , ■ ■ •, n n . 
To simplify the notation, we use mi, • • •, ton instead of n 1 , ■ ■ -,n n in what follows. The corresponding permutations 
are denoted by qi, ■ • •, q^. With these conventions, Eq. m becomes 


n\ ,ri2 ,••• ,11/v 

m T m 2)'") m lV 


J ( J 2 (- 1 ) Fpi ' P2 '-'^ n iIi Wn Pl {xi) 


\Pi,P2,--,Pn 


E (-1) 


(ii) Qx 1 <x2<---<x N e 2ikeoXj dx. (A2) 


\qi,q2,--,qN 


Since the integral in Eq. (IA2I) can be interpreted as a Fourier transform with respect to the coordinate Xj , we first 
evaluate the integral over the coordinates xi,- ■ ■ ■ ,Xn ■ Except for the sector function & Xl <x 2 <-<x N i 

the integrand in Eq. (IA2I) is symmetric under the exchange of any two variables that are smaller and larger than Xj. 
By changing the order of the integration variables that are smaller and larger than Xj , we get [53, |5|| [fjj, |7l| 


V 


i 

ni,n 2 , 
m i ,m 2 


,n N 

,m N 


P 1 ,P 2 ,---,PN 91.92 ,qN m qi ,m q2 ,-,m qN 


(A3) 


where 


T J . 

TYlq-^ , 771(22 


1 




dxj 


Vn P4 (Xj)Vm q4 (Xj)e 2lk ‘ oX i X 



(A4) 


with 


and 




( x)dx 



{x)dx. 


(A5) 


(A 6 ) 


According to the generalized Feldheim identity, the product of any number of Hermite polynomials can be expanded 
into a finite sum of Hermite polynomials |72 j, 


H Nl ( x) ■ ■ ■ H Nm ( x) 


E 


•••I'm 




1 H m (x), 


where 

771— 1 

M=J2 ( N i - 2 ui) + N m 
1=1 

and 


(A7) 


(AS) 


= n 


771—1 
1 = 1 


N, 


/+1 

n 


Ei=\ (N k -2u k + N t ) 
vi 


2 VI W- 


(A9) 
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The limits of the summation indices v k i n Eq. (1X71) are given by 

/ m -2 \ 

0 < vi < min (Ni,N 2 ) , 0 < v 2 < min (N 3 , Ni + N 2 - 2i/i) , —, 0 < v m -\ < min ( N m , 53 ( N k - 2 v k ) + N m _i j . 

(A10) 


k =1 


Using Eqs. (IA7I) - (IA10I) to rewrite the product of the two Hermite polynomials contained in I^\xj) and I^(xj), we 
obtain 


min (n Pk ,m qk ) 


' t'k . 1 y/c' p 

,m qk ( X j) = N n Pk Nr n qk 51 J_ 


and 


U fe) =o 


ain(n p ,m q ,) 


^„ Pfc+m 9 fe - 2 ,«( ;C ) eX P(- a: ) da! 


(All) 


v ' poo 

^ 2 p ) 1 ,m 9 ,fe)=^n Pi A m9! ^ <*„« / H _ 2 i/ m(x)exp(-x 2 )dx. 

n\ J Xn 


(A12) 


Instead of working with Eqs. (IA11I) and (IA12I) directly, we replace the Hermite polynomials in the integrals using the 
generating function g(x,t) of the Hermite polynomials, 


g(x,t) = e t2+2tx = V H n(x)'- 


n—0 


From Eq. (1A13I) . we have 


H n (x) = 


d n g(x,t) 


dt r ‘ 


(A13) 


(A14) 


t —o 


Inserting Eq. (1A14I) into Eqs. (1A1 1 1) - (1A121) . we find 


In p m q ( X j) = N n N m 2_^ a „( fe ) 7 (fc) 


r x i 

/ exp(—f 2 + 2 te — x 2 )dx 

J —oo 


(A15) 


i=0 


and 


In Pv m qi {Xj) — N npi Nm qi 51 


^ +m?i _ 2 „C0 T 


U !) =o 


9t n Pi+ m <!i- 2l/ l 


(0 


pOO 

/ exp(—t 2 + 2 te — x 2 )dz 

J Xn 


(A16) 


t=o 


The key point is that the Gaussian integrals can be calculated analytically. This yields 

min(n Pfe ,m Qk ) 


v[ k) =0 


J n Pk +m qk — 2u[ k) ,0 


UU„>>> = E »„<*> 

ffcl _ \ / 

+ exp(—Xj)H n ^ +m ^_ 2i> {ky_ 1 (xj)(l - 6. 


n Pk +m qk — 2v[ k) ,0 


(A17) 


and 


min(n Pi ,m qi ) 

In Pl ,m qi ( X j) = N n Pl Nm qi 5Z a i/<° 




^ n Pl +m qi —2v[ l \o 


exp( x j )H n ^ +m ^_ 2 i/ m_ 1 (xj)( 1 < 5 ripi+m?i _ 2 i ,( i >,o) 


, (A18) 
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where erf(a;) is the error function, 


2 f x 

erf(:r) = —= / exp(— t 2 )dt, 
v n Jo 


(A19) 


and S n , m the Kronecker delta function. Plugging Eqs. (IA17I) and (IA18I) into Eq. (IA4D and rearranging the order of 
the sums and products, Eq. (IA4I) can be written as a finite sum over one-dimensional integrals, 


nLiN Pk N mqk 


n Pi’ n P2’-’ n p N _ ivCAT — iV 


min(n P1 ,m qi ) 

E 


r 

min (n Pj _ 1 ,m g:j _ 1 ) 

f min (n p . +1 ,m, i+1 ) , 

i x ' 

•• E 

\ X E { 

l 

^- l5 =o 

{ 4 3+1) =o 1 


min {n PN ,m qN ) 

E 




: ^1 7 5 






a (j+i) x 


>, (A20) 


where 


7+ . .. ^O-i) v ^ +1 "> . .. jJ N h 

n P1 ,n P2 ,---,n PN 1^1 5 > ^1 > 5 1 


CJ /, # ( A ) 

? 

m 9l > m <22 

n, 


) — / #n Pj . {Xj)Hm qj {xj) exp(2ik so Xj Xj ) X 


A fc<j 

n z>j 


^ - ^erf(*i) ) <*, 


n Pfc +m, fe -2yJ ,0 


« n + exp(-^ 2 )7J ( fc) _ 1 (^-)(l - <5 

pfc yfc i 


n Pfe+ m 9fc _ 2 l/ l fc) " 0 


^ + ^erf(^))5 


. „ (i) n — exp(—. „ (!) ,(a;,')(l — <5 . „ (!) „) 

n Pl +m qi —2u\',0 J ' n p , +m 9i — 2u\' — 1 v n Pj +m 9i —2iq ' ,0' 


To simplify the notation, we define the index function d(k), 


d(k) = n Pk + m qk - 2 i/f } . 


dxj. (A21) 


(A22) 


For k < j and l > j, we use ki,k 2 , ■ ■ -,kic and h, I 2 , • • •, lc, respectively, to indicate all the k's and Vs that make d(k) 
and d(l) non-zero (0 < /C < j — 1 and 0 < £ < N — j). Then Eq. (IA21I) becomes 


r 


n P1 ’ n P2 i"'i n PN 

1^92 ) 




. (i-l) *,0+1) 

5 1 W \ 


,---A N) ) = 


J H np . ( Xj )H mq . ( Xj ) exp [2 iksoXj -{JC + C + l)x 2 ] (-1) £ - ^erf(^)^ + ^erf(^)^ 

H-d(ki) — 1 ip^j ) * * ’ -^d(fc;c )—1 )-^d(Zi) — 1 ) * * * ^-d{lc )—1 * (A23) 

/ , \ j—l — K, / , \ N—j — C 

Expanding the ( — -^erf(xj)J and ( + ^erf(xj)j terms, Eq. (|A23|) becomes 


N-j-C 


,b P\ i ,L P 2 > v fc PiV 
'WLq,TYlq2 ? '"'}77T.g^y 

N—C—1C—1 


0 - 1 ) ,, 0 + 1 ) 
i 


/ ii) _. 

K I 12 : 

j-X-K N-j-C 


. v Wx - 

1 "1 J ■ ■ Z/ 1 ) — 

(N-j-C 


(-i+ E E 


(— l) r x j exp [ 2 z/c so Xj — (/C + C + l)x^] [erf(xj)] r s x 

r=0 s=0 

J^ n P i ( x j)H mq i x j)Hd(ki)—l( x j) ‘ ‘ ' Hd(kic)-l( x j)Hd(l 1 )—l{Xj) • ‘ ‘ i (Xj)dXj. (A24) 


Equation (IA24D contains a product of Hermite polynomials, which can be converted into a finite sum of Hermite 
polynomials according to the Feldheim identity, 


JJn Pj { x j)H rnq ^{xj)H < i(j- 1 )_\(xj)---H ( i(j- K ^_x(xj)H ( i(i 1 )_i(xj)---H ( i{j c ^_x{xj) — ^ ^ o, l ^ 1 ... 1 /c+K+ 1 HM(xj ), (A25) 

V1---VC+K.+ 1 


where the coefficients a Vl ... vc+tc+1 and the relationship between M and the indices n Pi , rn qi , d(k\), d(lc) are defined 
in Eqs. (1X81)- (1A10I) . 
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Plugging Eq. (IA25I) into Eq. (1 A24I) . we find 


T l 


n P\ i n P 2 j” ‘ jv 

m n . 


K 


(i) 


. . . U- 1) jyO' + l) . . . 


(iV)\ 

) = 


tgi 5 , "'q2 ’ 

N—C—K—l 


j-l-ICN-j-C 


(-‘) £ E E 


r—0 s=0 


j-l-K\(N-j-C 


(-1 ) r E 


i'C+/C+i^(r+s;/C+£)> (A26) 


*'i---i'£+/c+i 


where 


r>M 

^(r+s;/C+£) 


/ oo 

exp [2zfc so a;j — (/C + £ + l)a: 2 ] [erf(^)] r+s H M (xj)dxj. 

-OO 


We call 0p:+ s -ic+c) the £ integral of the order r + s. The fy integral of order 0 can be evaluated analytically, 


■7 M 

A0,/C+£) 


V/c + £ +1 


exp - 


t 2 

v so 


/C + £ + l 


if 


2 ifc s 


M 


1 /(/C + £ + l)(/C + £) / V^ + >C + 1 


AC + £ 


M/2 


(A27) 


(A28) 


To evaluate £/(5 ? +s .jc+£) °f higher order, we develop an iterative procedure, in which the Q integral of order r + s is 
written in terms of Q integrals of order r + s — 1. Since the r + s = 0 integral is known, this allows for the evaluation 
of the Q integral of arbitrary order. Using the generating function of the Hermite polynomials, we have 


f qM roo 

G^+s-k+o = | frM J ex P [ 2( , ik so + t)xj - (1C + £ + l)a ; 2 - t 2 ] [erf(^)] r+s dxj j 
Using integration by parts, Eq. (IA29I) becomes 


t=o 


->M 


where 


and 


with 


r*i\ 

P( r+s-K+C) 


^(r+s,/c+£) — 'H{k S o,JC)£']t,Xj) [erf(xj)] 


( 5 M 

Jl) V 2 ) 

1 

|<9A m 

y(r+s,/C+£) y(r+s,/C+£) 

) 


t =o 


r+s |tCj —cxd 

Xj = — 00 


(A29) 


(A30) 


(A31) 


2 poo 

gfr+s,K+c) = ( r + s ) J K, £; t, Xj) [erf(a ; :) )] r+s_1 exp {-x 2 )dx 3 (A32) 


'H{k so ,1C,C-,t,Xj) = / exp [2(ik so + t)x — (AC + £ + l)a ; 2 — A 2 ] dx. 


The Gaussian integral in Eq. (1A33I) can be calculated analytically, 

k 2 so - 2 ik so t + {1C + C)t 2 


H{k ao ,K,,C\t,Xj) = 


2 V/C + £ + 1 


exp 


K. + C+ 1 


erf 


— ik so — t + (/C + £ -+■ 1) Xj 


V/c + £ + i 

(i) 


Using that the value of the error function erf(x) at infinity is known, erf(±oo) = ±1, g^ r+s /c+c) evaluates to 


,(!) 




g(r+sX+C) 2\/IC + £ + 1 
Inserting Eq. (IA34I) into Eq. (IA32I) . we have 

exp 


exp 


k 2 so — 2 ik so t + (AC + £)A 2 
/C + £ + 1 


[1 - (-l) r+a+1 ] . 


(A33) 


(A34) 


(A35) 


= (■■+») 1 ^ +£ ~ +1 


fc/„-2zfc so t+(A:+£)t 2 
/C+£+l 


erf 


— ik so — A T (/C + £ + l)x 


V/C + £ + 1 


(A36) 
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Next, we expand K+C) an( i 9{r+s K+c) terms of t up to power M and evaluate the Mth derivative with 
respect to t at t = 0. For /C + £ V 0, the exponential in Eq. (IA35f> can be interpreted as a generating function of the 
Hermite polynomials. Expanding this term into a sum of Hermite polynomials, we find 


o (1) 

y(r+s,/C+£) 


2 y/K + £ + 1 


[1 - (-l) r+s+1 ] exp 




/C + £ + 1 + £ + 1)(/C + £) V V /C + £ + l J \V /C + £ + 1 


/C + £ 


-t - 


K. + C 


2 V/C + £ + 1 


[1 _ (-ly+H-i] exp - 




/C + £ + 1 


E^« 

n=0 


2ifc s 




V(/C + £ + l)(/C + £) / n! 


For /C + £ = 0, the power series of g^+ s jc +C ) in terms of t is 

= [ 1 - ( - 1 ) r4 * tI ]”»(-*W 

This can be understood as the limiting result of Eq. (IA37I) , 


E 

n =0 


(2 iksoty 


n\ 


lim H„ 


2 ik„ 


K. + C 


n/2 


= (2 ik so ) n . 


. (A37) 


(A38) 


(A39) 


/c+£->o y -y/(/C + £ + 1 ) (/C + £) J + ^ + l 
In what follows, we take this limit when /C + £ = 0. The derivative of S , [r+s-A:+£) w ith respect to t at t = 0 is then 


9 M „(D \ 

Q^M ”(r+s;/C+£) J 


t=0 


2 V/C + £ + 1 



[• 2 

v so 


/C + £ + 1 




2 


M 


y(/C + £ + l)(/C + £) 



(A40) 


To evaluate f^ 2 ^ jc + £); we notice that it can be written in the form g^) +s ^ +/ r) / er f(''' )^ x j- The expansion of 
tc+C) §i ven in Eq. (1A37I) . An additional t-dependence enters through the error function in the integrand in 
Eq. (IA36I) . Rewriting the error function in a series in t, we have 


erf 


— ik so — t + (/C + £ + l)a 

VIC + £ + 1 


= erf ^ 


Vic + £ + ixn — 


ik. 


Vic + £ + i ) V^ \ic + £ +1 




AT / _ \ 

y exp [2ifc so a;j — (/C + £ + l):r 2 ] H n ( VIC + £ + lxj - 30 =j 

n=l ' V T T / 


1 


V/C + £ + 1 / \ VIC + £ + 1 / n! 


-. (A41) 


Using the multiplication theorem 


tf n (aaO = £ a”- 211 (a 2 - 1)“ (") ^-H n _ 2v {x) 

v=0 ' ' l ' 


(A42) 


and the addition theorem 


n 

H n (x + y) = ^2 

w =0 


H w (x)(2y) 


n—w 


(A43) 
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where [j] denotes the integer part of n/2, the Hernrite polynomial on the right hand side of Eq. (IA41I) can be expanded 
into a finite sum over products of Hermite polynomials in which the dependence on Xj has been “isolated”, 


H n ( V 1C + £ + 1 x ; — 


ik.« 


3 vic + c+i 

n—2v 


cr (l) M (-• -jy, (i) ^ w. (A44) 

D—0 w—0 X 7 \ / \ / 


For /C + £ = 0, Eq. (IA44I) contains an indeterminate term 0° which is understood to be 1. In this case, Eq. (IA44I) 
becomes 


H n (fij ik so ) — ^ ^ f j H w (2 ik so ) H n _ w (xj'). 

w=0 ' ' 


(A45) 


Using Eqs. (IA37I) . (IA41I) and (IA44I) . Eq. (IA36I) becomes 


a {2) 

VW+s-K+C) 


£■ 

n—0 


r + s 


VIC + C + 1 
2 (r + s) 


exp 




/C + £ + l 

n 111] u—2v 


Hn 


2 ikc: 


arEEE"”-* 


/C + £ 

V(/c+£ + i)(/c+£)y V^+^+l 

2 ik so \ .., n z u -\-v 


n/2 


L {r+s-l,K.+£) 


sfn(1C + £ + 1) 2 u=0 u=Q U1=0 


V(AC+ £ + !)(£ + £) 


(/C + £)~ (/C + £ + 1)“ 


n\ / u\ /v, — 2iA (2i>)! 
u J \ 2 v ) \ w J v\ 


-Hu 


2 ik s 


K. + C + 1 


so \ s-,(u—2v—w 


Q^Z-ix + c + i)\-r (A46) 


where 


I, 


(r+s—l,/C+£) 


/“ ertfy 

J — oo \ 


Vic + £ + iir 7 - — 




V/C + £ + 1 


[erf(a,’j)] r+s 1 exp (—x 2 )dxj. 


(A47) 


( 2 ) 

The Mth derivative of 5 ( r _ y a -ic+c) w ith respect to t at t = 0 is then 


d M ( 2 ) 

Q t M9(r+s-K+C) 


( 2 , 1 ) ( 2 , 2 ) 

S/r+s,AM-£) + #(r+s,/C+.C)’ 


(A48) 


t=o 


where 


o (2,1) 

»(r+s,/C+£) 


r + s 


V/C + £ + 1 


exp 




/C + £ + 1 


H 


2 ik s 


M 


JC + £ 


V(]C + C + l )(£ + £) / V^ + >c + i 


M/2 


"(r+s—l,/C+£) 


(A49) 


and 


( 2 , 2) 

y (r+s,t C+£) 


2 (r + s) 


M Ilf] u—2v 


V^{IC + £ + 1) M 2 + 


r £ £ £ 


2 ik s 


u —0 u—0 it/—0 


V , (A + £ + l)(K + £) 


u J \2vJ \ w J v\ 


(/C + C)—^ v (/C + C + l)” 3 ” x 

^ 8f;.-rxUi>- ( A5 °) 


/C + £ + 1 


To evaluate the integral I( r+S _!, k;+£), we integrate by parts. Using 

f [erf(x )]’ +s_1 exp(— x 2 )dx = 
Jo 


VZ[eri(x j )] r+ ‘ 
2 (r + s) 


(A51) 
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and 


rrffvr | r I It- lkso 

v^IerfO j)r +S 

Xj =00 /— 

l {' 1 1 J Vic + c + ij 

2 (r + s) 

x i= -oo 2 ( r + S ) 


[1 - (-l) r+5+1 ] , (A52) 


we find 


i (r+ -« +£) = ^ [i - ( -ir+-»] - 


v / A + £TTexp 


r + s 


J(r+s,K+C), (A53) 


where 


/ OO 

[erf(^)] r+s exp[-(/C + C + l)x| + 2 ik so Xj]dxj 

-OO 


(A54) 


Using Eq. (IA53I) in Eq. (1A49I) and using Eq. (IA40I) . we find 

nM _ st _ (2,2) 

b(,--M;A:+£) ~ J(r+s,K+C) ff(r+s,/C+£)' 


(A55) 


Our manipulations have reduced Q^+s-k+c) t° an expression that contains two types of one-dimensional integrals. 
The first one, J( r+St ic+C)i on ly depends on k so , 1C and C, implying that only a few of these integrals need to be 
calculated for a given k so . For r + s = 1, Jnx+C) can be calculated analytically, 


J 1CK+C) \j ^ + c + ! ex P ( /c + £ + i) erf 

For r + s > 1, J(\x+C) can be calculated numerically with essentially arbitrary accuracy. The second integral, 
G^+gZ^ic+c+i) ^ integral of order r + s — 1), is contained in x,+c)' see the iterative structure of our 
result more clearly, we rewrite Eq. (IA55j) as 


ik s 


^/(/C + £ + l)(/C + £ + 2) 


(A56) 


r>M _ sj \ ^ r>u—2v—w 

V{r+s-,lC+C) — U(r+s,/C+£) / , C J^(r+s-1-X.+C+l} ’ 

j 


(A57) 


where j runs over all combinations of allowed indices and Cj contains all the prefactors [the Cj can be read off 
Eq. (1A50I) ]. Thus, to determine the Q integral of order r + s, a finite number of Q integrals of order r + s— 1 is needed. 
To evaluate the Q integral of order r + s — 1, a finite number of Q integrals of order r + s — 2 is needed, and so on. 
Since the Q integral of order 0 is known analytically [see Eq. (IA28I) ]. the Q integral of arbitrary order can be obtained 
iteratively. It should be noted that the Q integral of order 1 has an analytical result since {?(o,it+£) and J{\ ,k+£) are 
known analytically. 

In summary, first we calculate, using Eq. (1A54I) . all the J integrals needed for the evaluation of the Q integrals. 
Secondly, using the values of the J integrals, the required Q integrals are calculated iteratively based on Eq. (1A55I) . 
Thirdly, plugging the values of the Q integrals into Eq. (IA26D , the T integrals with arguments ■ ■ ■, 1) , v ^ +1 ^, • • 

•, v[ N ^ are obtained. Forthly, plugging the values of the T integrals into Eq. (IA20I) . the I integrals can be calculated. 
Finally, plugging the values of the / integrals into Eq. (IA3I) . we obtain the values of the V integrals. The sums in the 
above expressions are all finite. As a result, the only error in evaluating the V> integrals comes from the numerical 
evaluation of the one-dimensional J integrals. Everything is analytical for N = 2 while one and three J integrals 
have to be evaluated numerically for N = 3 and 4, respectively. The reader may wonder why we choose the outlined 
iterative procedure for evaluating the one-dimensional integral given in Eq. (IA27I) over an evaluation of Eq. (IA27I) by 
direct numerical integration. The answer is two-fold. First, Eq. (IA27[) has to be evaluated for Hermite polynomials 
of different orders; the numerical integrals required in our iterative procedure, in contrast, are independent of the M 
quantum number. Second, for large M, the integrand in Eq. (IA27I) is highly oscillatory, making the one-dimensional 
integration somewhat non-trivial. 

To be concrete, we discuss how to apply the formalism to the three particle system. In this case, the V integral 
depends on j (j = 1,2,3) and the quantum numbers n\, n 2 , u. 3 , mi, m 2 , and m 3 . For a fixed j, there are 3! x 3! = 36 
terms (/ integrals) in the sum on the right hand side of Eq. (IA3I) . For each I integral, we have a double sum on 
the right hand side of Eq. (IA20I) . For j = 2, e.g., the summation dummies are and . Correspondingly, 
the T integrals have two arguments. The indices 1C and £ in Eq. (IA23I) satisfy 0 < YC < 1 and 0 < C < 1. As 
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a result, the indices r and s defined in Eq. (IA24I) satisfy 0 < r < 1 and 0 < s < 1 with the additional constraint 
r + s<N —1 — 1C — C, which is due to the limits r < j — JC — 1 and s < N — C — j in the sums on the right hand 
side of Eq. (IA24I) . Thus, the allowed (r + s,/C + C) combinations are (2, 0), (1, 0), (0, 0), (1,1), (0,1), (0, 2). Since the 
number of Hermite polynomials in the product in Eq. (IA24I) is JC + C + 2, which is less or equal to 4 (since 0 < K, < 1 
and 0 < C < 1), the coefficients on the right hand side of Eq. (IA25I) have at most three indices. 


Appendix B: Expressions for ( S x (x )) and {S z (x)) for N = 3 — 4 

For TV = 3, the expressions for C\ x , C 2x . and C z in Eqs. flT7l) - (H51) in terms of the coefficients C\ — C 4 in Eq. flCll) 


are 


c lx = (c 4 - C 3 )*(C 2 + C 4 ) + (Ci - c 3 yc 2 + c 4 y, 


(bi) 


c 2x = C\C 3 + C x Cl - \C 2 \ 2 - \C 4 \ 2 , 


and 


c z = i[(C 4 - C 3 )*(C 2 - C 4 ) - (Cl - C 3 )(C 2 - C 4 )*}. 
The expressions for ni x (x),n 2 x (x), and n z (x) in Eqs. ©-(USD are 


/ oo 

|-D(0, 1 , 2 )| 2 0 a:i<X 2 <a; 3 [( 5 (a; - xi) + S(x - x 3 )\dx, 

-CO 


(B2) 

(B3) 

(B4) 


and 


/ CO 

|D(0, 1 , 2)| 2 @ Xi<X2<X3 6(x - x 2 )dx, 

-OO 

/ OO 

|O(0,1, 2 )| 2 0 xi<a; 2 <a: 3 [( 5 (a; - x 4 ) - S(x - x 3 )]dx. 

-00 


(B5) 


(B 6 ) 


For TV = 4, the expressions for C 4x , C 2x , C 4z , and C 2z in Eqs. (I50l) - (l5l1) in terms of the coefficients C\ — C 43 in 
Eq. (l49l) are 

c lx = i [{C 2 + C 5 ) (Cl + Cs)* + (C 2 + C 5 y (C 4 + C a )\ 

+ [C 3 (c w + c 6 y + c; (C W + c 6 ) + c 4 (c 7 + c 9 y + c * 4 (c 7 + c 9 )] , (B 7 ) 

c 2x = i [(C 3 + c 4 ) (Cl + C 8 )* + (C 3 + C 4 )* (Cl + Cfc)] 

+ ±= [Cs (C 9 + C w y + C* 5 (C 9 + C 10 ) + c 2 (Cs + C 7 y + C * (Cs + C 7 )] , (B8) 


Ci x = - \(C 2 - Cs) (I Cl + C 8 ) 4 


-(C 2 -Cs)* (Ci + C s )] 

+ [Cs (C 10 - C 6 y - C* (C 10 - Cs) + c 4 (C 9 - C 7 y - C * 4 (C 9 - C 7 )\ , (B 9 ) 


and 


C 2z = - [(C7 3 - C 4 ) (Ci + C 8 ) 4 


-(Cs-cyf (Cr+Cg)] 

+ [Cs (c 9 - c w y - c 5 * (c 9 - Cw) + c 2 (c 7 - c 6 y - c * 2 (c 7 - c 6 )]. (bio) 
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The expressions for ni x (x),n 2 X {x),ni z (x), and 77 - 22 ( 2 ) in Eqs. (l50l) - (l5lT) are 

/ OO 

1-^(0? 1? 2,4)1 ©an<a;2<^3<^4 *^ 1 ) 4” 6(x ^ 4 )] dx, 

-oo 

/ oo 

|£>(0, 1, 2, 4)| 2 ©a?,<a; 2 <£c 3 <cc 4 [<*(x - X 2 ) + S(x ~ X 3 )]dx, 

-00 

/ oo 

\D(0, l,2,4:)\ 2 (d Xl<X2<X3<X4 [S(x - Xi) - 8{x - x 4 )\dx, 

-00 


and 


/ OO 

|-D(0, 1, 2, 4)1 ©xi<a;2<a;3<3;4 [^(*^ 22 ) d(2 X^^dx. 

-00 


(Bll) 


(B12) 


(B13) 


(B14) 
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